## Loading required package: assertthat
## Warning: package 'assertthat' was built under R version 4.0.3
##
## Attaching package: 'assertthat'
## The following object is masked from 'package:tibble':
##
## has_name
## Loading required package: glue
## Warning: package 'glue' was built under R version 4.0.3
##
## Attaching package: 'glue'
## The following object is masked from 'package:dplyr':
##
## collapse
## Loading required package: rlang
## Warning: package 'rlang' was built under R version 4.0.3
##
## Attaching package: 'rlang'
## The following object is masked from 'package:assertthat':
##
## has_name
## The following objects are masked from 'package:purrr':
##
## %@%, as_function, flatten, flatten_chr, flatten_dbl, flatten_int,
## flatten_lgl, flatten_raw, invoke, list_along, modify, prepend,
## splice
## `summarise()` has grouped output by 'ifbl_4by4', 'Jaar', 'TaxonIDParent'. You can override using the `.groups` argument.
##
## -- Column specification --------------------------------------------------------
## cols(
## id = col_double(),
## code = col_character(),
## naam_nederlands = col_character(),
## naam_wetenschappelijk = col_character(),
## taxn_vrs_key = col_character(),
## taxongroep = col_character(),
## usageKey = col_double(),
## scientificName = col_character(),
## rank = col_character(),
## order = col_character(),
## matchType = col_character(),
## phylum = col_character(),
## kingdom = col_character(),
## genus = col_character(),
## class = col_character(),
## confidence = col_double(),
## synonym = col_logical(),
## status = col_character(),
## family = col_character()
## )
## i Using ',' as decimal and '.' as grouping mark. Use `read_delim()` for more control.
##
## -- Column specification --------------------------------------------------------
## cols(
## substraat = col_character(),
## substrate = col_character()
## )
| phylum | n_obs | n_soorten | n_hokken |
|---|---|---|---|
| Anthocerotophyta | 218 | 4 | 123 |
| Bryophyta | 95624 | 410 | 796 |
| Marchantiophyta | 17486 | 109 | 727 |
## `summarise()` has grouped output by 'phylum'. You can override using the `.groups` argument.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Joining, by = "substraat"
| phylum | substrate | number of species |
|---|---|---|
| Bryophyta | artificial stone | 84 |
| Bryophyta | bryophyte | 3 |
| Bryophyta | decaying vegetation | 21 |
| Bryophyta | decorticated wood | 48 |
| Bryophyta | epiphytic on living wood | 73 |
| Bryophyta | floating on water | 1 |
| Bryophyta | gravel or sand | 108 |
| Bryophyta | mineral soil | 211 |
| Bryophyta | peat | 47 |
| Bryophyta | rock, hard | 95 |
| Bryophyta | rock, soft | 38 |
| Bryophyta | soil on rock | 68 |
| Marchantiophyta | artificial stone | 9 |
| Marchantiophyta | bryophyte | 17 |
| Marchantiophyta | decaying vegetation | 18 |
| Marchantiophyta | decorticated wood | 16 |
| Marchantiophyta | epiphytic on living wood | 15 |
| Marchantiophyta | floating on water | 3 |
| Marchantiophyta | gravel or sand | 34 |
| Marchantiophyta | mineral soil | 67 |
| Marchantiophyta | peat | 30 |
| Marchantiophyta | rock, hard | 14 |
| Marchantiophyta | rock, soft | 14 |
| Marchantiophyta | soil on rock | 20 |
| phylum | ell_l | number of species |
|---|---|---|
| Bryophyta | 1 | 9 |
| Bryophyta | 2 | 19 |
| Bryophyta | 3 | 37 |
| Bryophyta | 4 | 86 |
| Bryophyta | 5 | 125 |
| Bryophyta | 6 | 188 |
| Bryophyta | 7 | 189 |
| Bryophyta | 8 | 129 |
| Bryophyta | 9 | 15 |
| Marchantiophyta | 1 | 1 |
| Marchantiophyta | 3 | 20 |
| Marchantiophyta | 4 | 51 |
| Marchantiophyta | 5 | 51 |
| Marchantiophyta | 6 | 60 |
| Marchantiophyta | 7 | 56 |
| Marchantiophyta | 8 | 18 |
| phylum | ell_r | number of species |
|---|---|---|
| Bryophyta | 1 | 10 |
| Bryophyta | 2 | 63 |
| Bryophyta | 3 | 88 |
| Bryophyta | 4 | 67 |
| Bryophyta | 5 | 95 |
| Bryophyta | 6 | 141 |
| Bryophyta | 7 | 215 |
| Bryophyta | 8 | 101 |
| Bryophyta | 9 | 17 |
| Marchantiophyta | 1 | 28 |
| Marchantiophyta | 2 | 60 |
| Marchantiophyta | 3 | 20 |
| Marchantiophyta | 4 | 52 |
| Marchantiophyta | 5 | 30 |
| Marchantiophyta | 6 | 38 |
| Marchantiophyta | 7 | 21 |
| Marchantiophyta | 8 | 6 |
| Marchantiophyta | 9 | 2 |
| phylum | ell_f | number of species |
|---|---|---|
| Bryophyta | 1 | 12 |
| Bryophyta | 2 | 12 |
| Bryophyta | 3 | 61 |
| Bryophyta | 4 | 152 |
| Bryophyta | 5 | 187 |
| Bryophyta | 6 | 176 |
| Bryophyta | 7 | 51 |
| Bryophyta | 8 | 69 |
| Bryophyta | 9 | 54 |
| Bryophyta | 10 | 14 |
| Bryophyta | 11 | 4 |
| Bryophyta | 12 | 5 |
| Marchantiophyta | 4 | 15 |
| Marchantiophyta | 5 | 51 |
| Marchantiophyta | 6 | 52 |
| Marchantiophyta | 7 | 50 |
| Marchantiophyta | 8 | 47 |
| Marchantiophyta | 9 | 32 |
| Marchantiophyta | 10 | 9 |
| Marchantiophyta | 11 | 1 |
| phylum | ell_n | number of species |
|---|---|---|
| Bryophyta | 1 | 32 |
| Bryophyta | 2 | 195 |
| Bryophyta | 3 | 139 |
| Bryophyta | 4 | 196 |
| Bryophyta | 5 | 137 |
| Bryophyta | 6 | 71 |
| Bryophyta | 7 | 27 |
| Marchantiophyta | 1 | 43 |
| Marchantiophyta | 2 | 97 |
| Marchantiophyta | 3 | 45 |
| Marchantiophyta | 4 | 36 |
| Marchantiophyta | 5 | 23 |
| Marchantiophyta | 6 | 8 |
| Marchantiophyta | 7 | 5 |
| phylum | ell_t | number of species |
|---|---|---|
| Bryophyta | 2 | 52 |
| Bryophyta | 3 | 269 |
| Bryophyta | 4 | 142 |
| Bryophyta | 5 | 155 |
| Bryophyta | 6 | 108 |
| Bryophyta | 7 | 58 |
| Bryophyta | 8 | 13 |
| Marchantiophyta | 2 | 27 |
| Marchantiophyta | 3 | 113 |
| Marchantiophyta | 4 | 51 |
| Marchantiophyta | 5 | 43 |
| Marchantiophyta | 6 | 13 |
| Marchantiophyta | 7 | 2 |
| Marchantiophyta | 8 | 8 |
Only remove substrate levels with less than 5 species (for Ellenberg variables this is not necessary because they are included as continuous covariate).
## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
| phylum | substrate | n |
|---|---|---|
| Bryophyta | artificial stone | 84 |
| Bryophyta | decaying vegetation | 21 |
| Bryophyta | decorticated wood | 48 |
| Bryophyta | epiphytic on living wood | 73 |
| Bryophyta | gravel or sand | 108 |
| Bryophyta | mineral soil | 211 |
| Bryophyta | peat | 47 |
| Bryophyta | rock, hard | 95 |
| Bryophyta | rock, soft | 38 |
| Bryophyta | soil on rock | 68 |
| Marchantiophyta | artificial stone | 9 |
| Marchantiophyta | bryophyte | 17 |
| Marchantiophyta | decaying vegetation | 18 |
| Marchantiophyta | decorticated wood | 16 |
| Marchantiophyta | epiphytic on living wood | 15 |
| Marchantiophyta | gravel or sand | 34 |
| Marchantiophyta | mineral soil | 67 |
| Marchantiophyta | peat | 30 |
| Marchantiophyta | rock, hard | 14 |
| Marchantiophyta | rock, soft | 14 |
| Marchantiophyta | soil on rock | 20 |
## Warning: Ignoring unknown aesthetics: text
## `geom_smooth()` using formula 'y ~ x'
## Warning: Ignoring unknown aesthetics: text
## Warning in qlogis(prop_2000_2019): NaNs produced
## Warning in qlogis(prop_2000_2019): NaNs produced
## Warning in qlogis(prop_2000_2019): NaNs produced
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 7 rows containing non-finite values (stat_smooth).
## # A tibble: 126 x 2
## n_hok_1980_1999 n
## <int> <int>
## 1 0 245
## 2 5 34
## 3 6 18
## 4 7 18
## 5 8 12
## 6 9 22
## 7 10 12
## 8 11 14
## 9 12 11
## 10 13 8
## # ... with 116 more rows
## # A tibble: 142 x 2
## n_hok_2000_2019 n
## <int> <int>
## 1 0 18
## 2 1 87
## 3 2 60
## 4 3 35
## 5 4 32
## 6 5 26
## 7 6 27
## 8 7 12
## 9 8 17
## 10 9 25
## # ... with 132 more rows
To assess the changes in distribution of bryophyte species form 1980 until 2018 we divided the survey in two periods and compared the period 1980-1999 with the period 2000-2018. To minimize the effect of badly prospected grid cells we only used grid cells (16km²) where at least 40 species were found in both periods (Figure 2.1).
Figure 2.1: 16 km² grids of Flanders and the Brussels Capital Region with at least 40 species of Bryophytes recorded in both the period 1980-1999 and the period 2000-2019.
This resulted in 224 grid cells prospected in both periods. In these selected grid cells 35648 bryophyte observations (expressed as unique species per year and per 16 km² grid) were done in the first period and 45236 in the second period. To correct for potential changes in survey effort within those grids we used a relative change index as proposed by Telfer, Preston, and Rothery (2002), with slight modifications.
This method assumes that during the two surveys recorders attempted to observe as many species as possible in each grid cell. The size of the distribution area of individual species in a given period, i.e., the number of 16 km² grid cells where the species was observed, was expressed as a proportion of the total number of grid cells studied during that period. To summarise the changes in distribution area for all species, we modeled the number of occupied 16 km² squares in the second survey (out of the total number of 16 km² squares considered, i.e. proportional data) as a function of the logit-transformed proportion in the first survey (baseline) and covariates (phylum, substrate and Ellenberg values). Substrate effects were modeled for each phylum separately as random effects coming from a normal distribution with mean zero and an estimated standard deviation. Because the counts were overdispersed (variance larger than Binomial distribution), a beta-binomial distribution was used to account for overdispersion and a logit-link was used to link the proportions to the linear predictor with the following model equation:
\[ g(\mu_{i,s2}) = (\beta_0 + \beta_1\log\frac{\mu_{i,s1}}{1-\mu_{i,s1}} + \beta_L\text{ell}_L + \dots)*\text{phylum}_k + \mathbf{b}_{o,j} + \mathbf{b}_{1,j}\text{phylum}_k \] where:
\[ \begin{aligned} \mu_{i,s2} \equiv \mathbf{E}(Y_{i,s2}|\text{trials},\mathbf{b})&& \text{expected value}\\ g(\mu_{i,s2}) = \log\frac{\mu_{i,s2}}{1-\mu_{i,s2}} && \text{logit-link} \\ Y_{i,s2} \sim \text{BetaBinomial}(\mu_{i,s2}, \phi_{i,s2}) && \text{Beta-Binomial distribution}\\ \begin{bmatrix}\mathbf{b}_{o,j}\\ \mathbf{b}_{1,j}\end{bmatrix}\sim N(0,\Sigma) && \text{random intercept and slope}\\ \Sigma = \begin{bmatrix}\sigma^2_{0,j}&0\\0&\sigma^2_{1,j}\end{bmatrix} && \text{variance-covariance matrix} \end{aligned} \]
Species that occurred in less than five 16 km² grid cells in the first period were omitted form the analyses. We also used only species that were present in both periods.
A Bayesian hierarchical modeling procedure was used to calculate the posterior distribution for these parameters (Bürkner 2017) and we assumed non-informative priors.
We focused model selection on the Ellenberg terms only. All other effects were kept in the model. First, interaction terms for Ellenberg-values with phylum for which the 95% confidence interval contained 0 were removed. Next, main Ellenberg terms for which the 95% confidence interval contained 0 were removed. The simplest model with the highest predictive accuracy based on leave-one-out cross validation (Vehtari, Gelman, and Gabry 2017) was selected.
After rearranging the model equation (for one phylum) and exponentiating the linear predictor, ignoring all random effects, we obtain the following odds-ratio:
\[ \frac{Odds_{\mu_{i,s2}}}{Odds^{\beta_1}_{\mu_{i,s1}}} = \exp(\beta_0 + \beta_L\text{ell}_L + \dots) \] which allows us to predict the effect of Ellenberg values on the relative increase or decrease of the Bryophytes.
In a similar way, substrate random effects are calculated from the posterior linear predictor, conditional on Ellenberg values fixed at their median values. Gelman, Hill, and Yajima (2012) point out that no correction for multiple comparisons are necessary when substrate effects are calculated in this way.
All analyses were done using the R language for statistical computing (R Core Team 2020).
## Start sampling
## Start sampling
## Warning: Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.
| elpd_diff | se_diff | elpd_loo | se_elpd_loo | p_loo | se_p_loo | looic | se_looic | |
|---|---|---|---|---|---|---|---|---|
| alt_telfer_3 | 0.00 | 0.00 | -3316.11 | 31.41 | 22.11 | 3.34 | 6632.21 | 62.82 |
| alt_telfer_2 | -0.88 | 1.27 | -3316.99 | 31.10 | 23.61 | 3.27 | 6633.98 | 62.20 |
| alt_telfer | -2.67 | 1.52 | -3318.78 | 31.16 | 26.16 | 3.75 | 6637.55 | 62.32 |
plot(loo1$loos$alt_telfer, label_points = TRUE)
plot(loo1$loos$alt_telfer_3, label_points = TRUE)
alt_telfer_3$data[745,]
## n_hok_2000_2019 n_hok_tot_2000_2019
## 745 8 228
## qlogis(n_hok_1980_1999/n_hok_tot_1980_1999) phylum ell_f ell_t
## 745 -1.029619 Marchantiophyta 4 3
## n_hok_1980_1999 n_hok_tot_1980_1999 substrate
## 745 60 228 epiphytic on living wood
analyse_data_telfer_agg[745,]
## # A tibble: 1 x 14
## phylum taxon_id_parent n_hok_1980_1999 n_hok_2000_2019 n_hok_tot_1980_1~
## <chr> <dbl> <int> <int> <int>
## 1 Marchantiop~ 760 60 8 228
## # ... with 9 more variables: n_hok_tot_2000_2019 <int>, prop_baseline <dbl>,
## # prop_2000_2019 <dbl>, ell_l <dbl>, ell_n <dbl>, ell_f <dbl>, ell_r <dbl>,
## # ell_t <dbl>, substrate <chr>
taxa[taxa$id == 760, ]
## # A tibble: 1 x 15
## id code naam_nederlands naam_wetenschap~ naam_duits naam_engels naam_frans
## <int> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 760 ptil~ Boomfranjemos Ptilidium pulch~ <NA> <NA> <NA>
## # ... with 8 more variables: creation_date <dttm>, creation_user <chr>,
## # update_date <dttm>, update_user <chr>, taxn_vrs_key <chr>,
## # taxon_groep_id <int>, naam <chr>, omschrijving <chr>
#Boomfranjemos looks like an outlying observation
#it is the only epiphytic liverworth on living wood that declines strongly:
analyse_data_telfer_agg %>%
filter(substrate == "epiphytic on living wood",
phylum == "Marchantiophyta")
## # A tibble: 8 x 14
## phylum taxon_id_parent n_hok_1980_1999 n_hok_2000_2019 n_hok_tot_1980_1~
## <chr> <dbl> <int> <int> <int>
## 1 Marchantiop~ 450 142 213 228
## 2 Marchantiop~ 558 170 146 228
## 3 Marchantiop~ 559 219 218 228
## 4 Marchantiop~ 591 93 204 228
## 5 Marchantiop~ 600 5 13 228
## 6 Marchantiop~ 760 60 8 228
## 7 Marchantiop~ 778 99 193 228
## 8 Marchantiop~ 842 14 4 228
## # ... with 9 more variables: n_hok_tot_2000_2019 <int>, prop_baseline <dbl>,
## # prop_2000_2019 <dbl>, ell_l <dbl>, ell_n <dbl>, ell_f <dbl>, ell_r <dbl>,
## # ell_t <dbl>, substrate <chr>
## Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
## TMB was built with Matrix version 1.3.2
## Current Matrix version is 1.2.18
## Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package
## Registered S3 method overwritten by 'broom.mixed':
## method from
## tidy.gamlss broom
## Warning in tidy.brmsfit(alt_telfer_3): some parameter names contain underscores:
## term naming may be unreliable!
| effect | component | group | term | estimate | std.error | conf.low | conf.high |
|---|---|---|---|---|---|---|---|
| fixed | cond | NA | (Intercept) | 0.4904610 | 0.1315195 | 0.2316668 | 0.7477269 |
| fixed | cond | NA | qlogisn_hok_1980_1999Dn_hok_tot_1980_1999 | 0.8964563 | 0.0173138 | 0.8631410 | 0.9301981 |
| fixed | cond | NA | phylumMarchantiophyta | -0.7365657 | 0.1734234 | -1.0744809 | -0.3895969 |
| fixed | cond | NA | ell_f | -0.0594441 | 0.0141049 | -0.0868474 | -0.0316259 |
| fixed | cond | NA | ell_t | 0.0265169 | 0.0186195 | -0.0099043 | 0.0630969 |
| fixed | cond | NA | qlogisn_hok_1980_1999Dn_hok_tot_1980_1999:phylumMarchantiophyta | 0.0336040 | 0.0396186 | -0.0438984 | 0.1119242 |
| fixed | cond | NA | phylumMarchantiophyta:ell_t | 0.0925428 | 0.0381588 | 0.0181477 | 0.1674919 |
| ran_pars | cond | substrate | sd__(Intercept) | 0.1613464 | 0.0518312 | 0.0851389 | 0.2810960 |
| ran_pars | cond | substrate | sd__phylumMarchantiophyta | 0.1092055 | 0.0826154 | 0.0042748 | 0.3041457 |
## Adding missing grouping variables: `phylum`
## Note: Using an external vector in selections is ambiguous.
## i Use `all_of(ellenberg)` instead of `ellenberg` to silence this message.
## i See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.
## Adding missing grouping variables: `phylum`
## Adding missing grouping variables: `phylum`
## Adding missing grouping variables: `phylum`
Bürkner, Paul-Christian. 2017. “Brms : An R Package for Bayesian Multilevel Models Using Stan.” Journal of Statistical Software 80 (1): 1–28. https://doi.org/10.18637/jss.v080.i01.
Gelman, Andrew, Jennifer Hill, and Masanao Yajima. 2012. “Why We (Usually) Don’t Have to Worry About Multiple Comparisons.” Journal of Research on Educational Effectiveness 5 (2): 189–211. https://doi.org/10.1080/19345747.2011.618213.
R Core Team. 2020. R: A Language and Environment for Statistical Computing. Manual. Vienna, Austria: R Foundation for Statistical Computing.
Telfer, Mark G, C. D. Preston, and Peter Rothery. 2002. “A General Method for Measuring Relative Change in Range Size from Biological Atlas Data.” Biological Conservation 107 (1): 99–109. https://doi.org/10.1016/S0006-3207(02)00050-2.
Vehtari, Aki, Andrew Gelman, and Jonah Gabry. 2017. “Practical Bayesian Model Evaluation Using Leave-One-Out Cross-Validation and WAIC.” Statistics and Computing 27 (5): 1413–32. https://doi.org/10.1007/s11222-016-9696-4.